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Abstract. We study quantum dynamics of bosonic atoms that are excited to form 
a phase kink, or several kinks, by an imprinting potential in a one-dimensional 
^ , trap. We calculate dissipation due to quantum and thermal fluctuations in soliton 

^ ' trajectories, collisions and the core structure. Single-shot runs show weak filling 

of a soliton core, typically deeper solitons in the case of stronger fluctuations and 
| spreading/disappearing solitons due to collisions. We also analyze a soliton system in 

an optical lattice that shows especially strong fluctuation-induced phenomena. 
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Improved technology to manipulate ultra-cold atomic gases by optical lattice potentials 
PQ, atom chips [2], magnetic fields, etc., has provided tools to trap atoms in strongly 
confined geometries and to adjust the dimensionality of the system. Such many-particle 
systems only interact weakly with their environment and, over short time-scales, they 
may in many situations be approximately considered as isolated and studied by models 
based on unitary quantum evolution. Moreover, the many-body atomic states may 
be engineered and controlled at high accuracy and the detailed dynamics of atoms in 
response to external perturbations can be experimentally investigated. 
^ ■ In a tightly-confined one-dimensional (ID) tube-potential quantum fluctuations 

of bosonic atoms can become strong [3], HI [5]. The quantum effects may be further 
enhanced, e.g., by applying an optical lattice potential along the axial direction of the 
trap [6j[7J|8]. Stochastic phase-space methods provide a useful description of dissipative 
bosonic atom dynamics in ID traps due to quantum and thermal fluctuations of the 
atoms when the particle number is not too low [9J, UTil HH [12] . The approximate quantum 
dynamics of atoms can be modelled in the truncated Wigner approximation (TWA) 
[131 [HI [15l HE] by unitary quantum evolution where the quantum statistical correlations 
of the initial state may be accurately synthesized for a variety of quantum states in 
the Wigner representation. The method can incorporate a very large phase-space, with 
a large number of degrees of freedom, in which case the atom dynamics, e.g., does 
not need to be restricted to the lowest energy band of an optical lattice. Dissipative 
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dynamics emerges from a microscopic treatment of the unitary quantum evolution 
without any explicit damping terms and there is no need for the frequently problematic 
separation of quantum dynamics to 'system' degrees of freedom and 'environment' [17J. 
In addition, it is possible to study excitations of the system far from thermal equilibrium, 
quantum expectation values of various experimental observables and outcomes of single- 
shot measurements. Previous unitary TWA simulations [10] were qualitatively able to 
produce the experimentally observed damping rate of a dipolar centre-of-mass motion 
of bosonic atoms in a very shallow, strongly confined ID optical lattice \7\. It was 
found that, due to quantum fluctuation induced momentum uncertainty, a small atom 
population reaches a critical velocity, representing an onset of a dynamical instability 
in the corresponding classical system. The calculated damping rate was approximately 
proportional to the atom population in the dynamically unstable velocity region. 

Here we study dissipative quantum dynamics and relaxation of a ID bosonic atomic 
gas that has been excited to form a moving phase kink by an optical imprinting potential. 
Phase kinks have been experimentally imprinted in this way in 3D atomic Bose-Einstein 
condensates (BECs) by imaging the atom cloud, e.g., through an absorption plate 
[T8| IT9] or by using optical holograms |20j. We recently showed how such a set- 
up could be used to study dynamics of dark solitons in situations where quantum 
fluctuations are important [21]. The integrability of the soliton dynamics is broken 
by a trap and quantum and thermal fluctuations promote sound wave emission that 
may dissipate and eventually equilibrate with the soliton. Numerically tracking soliton 
coordinates in individual stochastic realizations provided us with a tool to calculate 
quantum mechanical expectation values and uncertainties of the soliton position. We 
found, e.g., a surprising result that the quantum expectation value of the speed of a 
soliton is reduced due to enhanced quantum fluctuations, as a result of the nonlinear 
dependence of the soliton speed upon its phase distribution. Single-shot runs in an 
optical lattice revealed effects of dynamical instabilities, such as jittering oscillatory 
motion, splitting and disappearing solitons. In this paper, we study a detailed structure 
of a quantum soliton, its quantum statistics and dissipative dynamics due to quantum 
and thermal fluctuations in fairly large atom-number systems. In particular, we consider 
the effects of dissipative dynamics on a soliton core structure, predicting a distinct 
difference between the filling behaviour of a soliton core in single-shot experimental runs 
and in the quantum expectation value of atom density found by averaging over many 
experimental runs. Single-shot realizations show only weak filling of a soliton core at 
later times and typically deeper solitons for the case of stronger quantum fluctuations. 
On the other hand, quantum expectation values for the atom density and pair-correlation 
functions are smeared out due fluctuation-induced drifting of solitons along the trap - 
this is particularly visible for very slow solitons in an optical lattice system we consider. 
We also simulate quantum effects of collisions of soliton pairs and relaxation in a system 
of several solitons that may behave chaotically in the corresponding classical case. As 
in [2T] quantum and thermal fluctuations are synthesized within the truncated Wigner 
approximation in the quasi-condensate description. Dark solitons in nonlinear optical 
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fibers were previously modelled using TWA [22J. In atomic BECs, finite temperature 
dissipation has now been theoretically shown to manifest itself in an increasing oscillation 
amplitude of solitons in a harmonic trap [231 EH ESI EI] • Quantum properties of dark 
solitons [TFJ have recently started attracting considerable interest, including studies, 
e.g., of soliton core structure [261 EH EH1 EHl USB EH [32], phase kinks in a uniform space 
close to the Mott transition [33], and soliton statistics [33]. 

2. Dark soliton experiments 

Dark solitons have been prepared in many experiments, and their subsequent dynamics 
probed. The first dark solitons in ultra-cold atomic Bose gases were generated by 
imprinting a sharp phase jump across the gas by optical potentials [HI [T9] . or by 
imprinting density defects using slow light pulses [35l [36] . Subsequently, dark solitons 
have been created, e.g., by merging of two BECs [37] and by superfluid flow around a 
barrier [38]. In typical experiments the atoms were confined in 3D traps, or in ID traps 
with only a weak radial trapping frequency, and the systems are accurately described 
by the classical Gross-Pitaevskii equation (GPE). In the phase imprinting method, a 
soliton is generated by applying a constant 'light-sheet potential', of value V^, to a 
part of the atom cloud, for time r [TBI EH EO] • The light-sheet potential is obtained 
by shining a far-detuned laser, e.g., through an absorption plate, so that the resulting 
dipole potential for the atoms exhibits a sharp edge. In the classical GPE limit, the 
potential imprints a phase jump 



generating a dark (or grey) soliton (phase kink) with velocity v and the density dip at 
the phase kink (soliton core) with the minimum density p s [39] 



where c is the speed of sound and p& is the density of the atomic background. The 
soliton is stationary (dark) for <p c = ir, with a zero density at the kink. Other phase 
jumps produce moving (grey) solitons, with non- vanishing densities at the soliton core, 
so that |v| — > c for C — y 0. Soliton imprinting is accompanied by a density (sound) 
wave moving in the opposite direction to the soliton with speed approximately equal to c 
[T8] . created as a by-product of the perturbation of density by the imprinting potential. 

More complex defect structures may generally be imprinted on superfluids using 
optical phase holograms to shape laser fields, so that, via coupling with matter waves, 
the light acts as a hologram to shape the BEC [30]. A spatial light modulator was 
experimentally used to generate a multi-step potential and to prepare multiple phase 
kinks [JT]. The potential 



4> c = V^r/h, 



(1) 



v/c = cos(0 c /2), 
p s = p b cos 2 (4> c / 2) 



(2) 
(3) 



Vi(x, t) = Yl y^fa - t)6(x - Xj ) 



(4) 



j 
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(where 9 is the Heaviside step-function) imprints phases 4>j at position Xj. The potential 
Vfa may be negative, such that V^. = —V^, imprints a phase of 2tt — (f>j, corresponding 
to a soliton travelling in the direction of the negative x-axis. Hence the initial positions 
and velocities of solitons may be accurately controlled. Other experiments have created 
multiple solitons by colliding two BECs [121 S3] ■ 

Subsequent to their formation, single solitons in a harmonic trap oscillate at angular 
frequency to / a/2, where to is the frequency of the trap [M] • Although many soliton 
experiments were performed in harmonic traps, these had a 3D character [45] . such 
that solitons could decay into a hybrid of vortex lines and rings ('snake' instability) 
[46J; consequently, lifetimes of solitons were limited to less than one oscillation period. 
However, some experiments suppressed the snake instability by increasing the radial 
trapping potential [20l [43] , permitting observations of oscillations at the predicted 
frequency [20]. Recently, experiments probed the interactions between dark solitons 
[20[ W\\ 143] or between dark and dark-bright solitons [20] . Collisions between two dark 
solitons are accompanied by position-shifts [4T1 14"8] . which change the soliton oscillation 
frequency [13]. During experiments, solitons repeatedly collided up to seven times 
without losing their integrity [43J. 

Although some soliton experiments were sufficiently ID to suppress snake 
instability, residual radial dynamics affected the oscillation frequencies by approximately 
3% [43J. In more tightly confined ID traps the radial density fluctuations may be 
completely suppressed. For example, in an atom transport experiment, a 2D optical 
lattice divided a BEC into an array of decoupled ID tubes [7] with radial and axial trap 
frequencies in each tube u r = 2tt x 38kHz and u r = 2tt x 60Hz, respectively, and about 70 
atoms in the central tube. A wide variety of ID trapping schemes are also provided on 
atom chips using electromagnetic fields [2]. In a tightly confined ID limit, bosonic atoms 
become more strongly interacting at low atom density [19] , characterized by the effective 
interaction parameter 7i nt = mg/h 2 p, where p is the ID atom density. For increasing 
values of 7 int , enhanced phase fluctuations destroy the long-range phase coherence across 
the quasi-condensate, but density fluctuations remain suppressed-eventually the atoms 
reach the fermionized Tonks-Girardeau regime [61 H] . 



3. Truncated Wigner approximation 

In our analysis, quantum and thermal fluctuations of atoms are approximately included 
within TWA [131 El [EJ EEU HE]. Our TWA formalism and noise generation closely 
follows [11 i], except that here we fix the total atom number and use a quasi-condensate 
description for quantum statistical atom correlations |21j. We replace the quantum field 
operators (ip,^) by the classical fields (V'W) V'wOj governed by a ID unitary evolution 

d h 2 d 2 

ih—*l) W = {- ^Q^ + V + gNtot^wl^w, (5) 

where the interaction strength g = 2hu r a, the s-wave scattering length a, the total 
number of atoms At ot , and V = V ext + V% is the external ID potential including the 
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time-dependent imprinting potential. Equation (J5J) formally coincides with the classical 
GPE, but here quantum and thermal fluctuations are sampled in the initial state by an 
ensemble of stochastic fields ipw{x, 0) in the Wigner representation. Initially, before the 
phase imprinting, the atoms are assumed to be in thermal equilibrium in a harmonic 
trap, or in a combined harmonic trap and optical lattice. In a tight ID trap we expect 
significant phase fluctuations that we found particularly important to phase kinks [21] 
and in order to incorporate these we introduce a quasi-condensate description for the 
field operator [50] 



i>(x, 0) = \Jpo{x) + 6p(x) exp(i<p(x)). (6) 

The density 8p(x) and phase <p(x) operators may be written in the Bogoliubov-type 
expansion [50], requiring Sp/po and |<5ZA(/3| to be much less than one, where 81 is the 
spacing on the numerical grid on which we calculate the operators and A0 is the gradient 
of the phase operator across one gridpoint. Thus (for j > 0) 

^ = r~~rr H ( ( Pi( x ) & j ~ V*j( x )°i) > ( 7 ) 

1 ipo{x) j 



5p{x) = ypoOr) ( 6 Pj( x )®j + , (8) 

3 

where ipj(x) = Uj(x) + Vj(x) and Spj(x) = Uj(x) — Vj(x) are given in terms of the 
solutions to the Bogoliubov equations 

f-h 2 d 2 



+ V-n + 2N g \tfj \ 2 j Uj - Nogipfa = e jUj , (9) 
VV - n + 2iV#o| 2 ) v d - N gr 2 u 3 = -e jVj , (10) 



y 2m dx 2 

f-h 2 d 2 

\ 2m dx 2 

where p, is the chemical potential. Here po = iVolV'o^)! 2 and ipo(x) is the ground state 
wavefunction with Nq particles. 

We fix the total atom number to N tot and allow the ground-state and excited- 
state populations, No and iV nc , to fluctuate. In the Bogoliubov approximation, the 
expectation value of the excited state population (j > 0) 



N nc = JdxJ2 (\u 3 \x)\ 2 + |^(x)| 2 ) + |^-(x)p 



where (djotj) = flj = [exp(ej/fcsT) — l] -1 . In order to sample the stochastic initial 
state of TWA with the correct quantum statistics, the operators (aj, a>j) in (jTJ) and (jH]) 
are replaced by complex variables (a*,ctj), which may be generated from the Wigner 
distribution for harmonic oscillators with energy ej at temperature T [H] , resulting in a 
stochastic Wigner representation (<pw(x),5pw(x)) of phase and density operators. The 
stochastic initial state for the time evolution flSD then reads 



ipw{x, 0) = \Jpo,w(x) + Spw(x) exp(i(p w (x)), (12) 

where p ,w w ih be specified shortly. Due to the Wigner representation, the stochastic 
variables correspond to symmetrically ordered expectation values of operators (dj, acj) 
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and, e.g., the expectation values of the mode occupations differ from those calculated 
in the Wigner representation by half an atom per mode: = (a*aj)w — 1/2, 

where the notation {...)w indicates the expectation values are calculated in the Wigner 
representation. The excited state population for each run therefore reads 

iv nc = JdxY, [(hf - 1/2) (M*)| 2 + M*)| 2 ) + M*)| 2 ] , (is) 
3 

fluctuating around the correct expectation value (ITT]) . We set the ground-state 
population for each stochastic realization as No = N tot — iV nc , so we can set for the 
initial state in each run p 0jW (x) = (N + l/2)|-0 o (a;)| 2 . Then N fluctuates around 
Nq = Ntot ~~ N nc as required. 

In simulations we vary the ground-state depletion N nc /N tot . At temperature T = 
we keep the nonlinearity fixed at N tot g = lOOhul, but adjust the ratio g/N tot where 
I = (h/mu) 1 / 2 . This is tantamount to varying the effective interaction strength 7i n . We 
also study the effects of thermal depletion by varying T, by fixing both N tot and g and 
setting g/Nt t to be sufficiently small so that quantum fluctuations are not dominant. 

After evolution of the stochastic fields ipw, we extract physical quantities from 
the ensemble by calculating their quantum mechanical expectation values. Due to 
the symmetric ordering of the expectation values in the simulation data, we need to 
transform these to the normally ordered expectation values, corresponding to typical 
physical measurements [9J [UJ. We find that the soliton position measurements 
considered in this paper are unchanged in the normal-ordering of phonon-mode 
occupation numbers, but that densities, pair-correlations and number fluctuations are 
quantitatively affected by the ordering. Hence the soliton position coordinates in 
individual experimental runs can be accurately extracted from the Wigner density |"0w| 2 - 



4. Isolated soliton in a harmonic trap 

4-1. Damping and dissipation in oscillatory dynamics 

We study the imprinting of single phase kinks and their subsequent dynamics in a 
harmonic trap in TWA. In our previous study [21] we showed how quantum expectation 
values and uncertainties could be calculated for the dynamical variables of the soliton. 
We ran ensembles of hundreds of realizations (e.g. 400), and numerically tracked 
the soliton position coordinates x(t) from l^wl 2 in each run. This allowed us to 
calculate, e.g., the quantum mechanical position expectation value (x) and uncertainty 
5x = J (x 2 ) — (x) 2 . Based on our findings, we made, e.g., the following observations that 
are relevant to our present study: (i) soliton trajectories are damped by both quantum 
and thermal fluctuations, increasing the amplitude of soliton oscillation; (ii) quantum 
mechanical soliton position uncertainties increased as a function of time for systems 
with large quantum depletion A^ nc /iV tot -this increase was due to the initial uncertainty 
in soliton velocity and deviations in the soliton trajectories as they interacted with sound 
waves; (iii) one of our most dramatic findings was that the quantum expectation values 
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Figure 1. (a) Quantum expectation value of the damping parameter 7 of soliton 
oscillation amplitude in a harmonic trap vs temperature for fixed giVtot = 100^' an d 
N to t = 900. iV nc /iVtot varies between 0.01 at T = 0, and 0.2 at T = 22huj/k B ; (b) 
7 vs depleted fraction N nc /Ntot at T = for fixed gNtot when we vary the effective 
interaction g/N% Q t- No points are plotted after N n c/Ntot = 0.01 corresponding to 
Mot = 900 in which case fitting of the damped oscillations became less accurate. The 
Bogoliubov limit N nc /Ntot corresponds to 7/w ~ —0.023 that is different from 
the classical GPE result 7/cj ~ —0.014 (shown as an open red circle). The damping 
increases (7 becomes more negative) faster than linearly with depletion. 



for the speed of a soliton were lower in systems with large quantum depletion which 
is attributable to a broad distribution of phase jump values across the kink caused by 
quantum fluctuations. We mapped the phase kink of each trajectory in the ensemble to 
a corresponding speed using the classical formulae (J3]) and found that due to the negative 
curvature of the soliton speed | cos(0/2)|, a symmetric phase distribution always leads 
to a lower quantum expectation value for speed (| cos(0/2)|) than the speed given by the 
phase expectation value | cos((0/2))|. Consequently, a greater phase uncertainty leads 
to a lower expectation value for the soliton speed. 

Here we evaluate in more detail the dissipation in soliton dynamics due to quantum 
and thermal fluctuations. We also study the structure of the soliton core and phase in 
individual stochastic realizations of TWA that represent possible outcomes of single- 
shot experiments. We then compare the quantum statistics of a soliton core in single- 
shot outcomes to quantum expectation values of the atom density which would be 
obtained in an experiment by averaging the density over many runs. We simulate the 
optical imprinting of phase kinks according to ([5]) in the harmonic trapping potential 
Kx P = mcu 2 x 2 /2, with the imprinting potential @, for = 4166.7hu, applied over a 
period of r = 4.78 x 10~ 4 /u. Such a potential imprints a phase kink of <j) c = 2 in the 
classical (GPE) limit. Due to the stochastic sampling of quantum and thermal noise in 
the initial state of each realization of the Wigner wavefunction, the actual phase jump 
values across the soliton core fluctuate from run to run, resulting in a variation in the 
shape of the solitons and in their initial velocities. After the imprinting the solitons 
oscillate in the harmonic trap. 

Previously [2T], we observed that soliton oscillations in a harmonic trap exhibit 
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Figure 2. Individual stochastic realizations of single soliton dynamics in a harmonic 
trap at T = for fixed nonlinearity gNtot = lOOTuol. In (a)-(c) the Wigner density 
\ipw(x, t)\ 2 is plotted with iVtot = 50 and the imprinted phase in the corresponding 
classical case (j> c — 2.6 (a) and 4> c — 2.0 [(b)-(c)]. (d)-(f) correspond to the same 
stochastic trajectories as (a)-(c) [the same realizations of (a*, ctj)] but with iVtot = 100. 
Soliton interactions with sound waves increase dissipation and oscillation amplitude. 
Dissipation is enhanced by strong quantum fluctuations [compare (a)-(b) to (d)-(e)]. 



damping; since dark solitons have negative mass and kinetic energy, damping of a 
soliton increases its amplitude of oscillation. We track each soliton's trajectory by 
locating the local density minimum around the phase kink in ipw for each realization. 
We fit individual trajectories with the curve x(t) = f(t) exp(— jt), where f(t) is a 
sinusoid and 7 < 0. We separately study the damping due to the ground-state depletion 
resulting from quantum fluctuations (at T = 0) and from finite-temperature atoms (for 
parameters for which the corresponding T = damping is weak). We only fit trajectories 
for parameters for which fluctuations are not too large, such that the trajectories are 
sufficiently close to sinusoidal. Figures [1] (a) and (b) show that there is a faster-than- 
linear increase in damping with depletion in both cases. 

The source of this damping behaviour can be deduced by consideration of individual 
trajectories in ensembles with large depletion (shown in figures [2] and |3]). Strongly 
dissipative behaviour occurs when energy is exchanged between the solitons and sound 
waves. This is particularly clear for large quantum depletion shown in figures E^a)-(c). 
We find that at finite temperature large thermal depletion has an associated sound 
wave background [figure [3fd)- (f )] , but the soliton trajectories are less perturbed by the 
sound waves than in cases of large quantum depletion. We contend that the breaking of 
integrability by the harmonic trap causes the energy to disperse amongst the excitations. 
Solitons with relatively less negative kinetic energy are more prone to dissipation, and 
the damping rate is greatest in these trajectories. 

The limit g/N tot — > 0, for fixed gN tot , corresponds to N nc /N tot — > and we expect 
the Bogoliubov approximation to become increasingly accurate (if the ground state 
and the excited state populations are not solved self-consistently, N nc is approximately 
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Figure 3. Individual realizations of trajectories for a single imprinted soliton at 
different temperatures for iVtot = 900, gNtot — WOHujI and C = 2.0. The Wigner 
density \^ w {x, t)\ 2 is plotted for (a) T = 0, (b)-(c) T = 16hu)/k B , and (d)-(f) 
T = TlhbjjkB- For large thermal depletion, sound wave emission is clearly visible, 
as in the case of large quantum depiction in figure [5] Interactions between the sound 
waves and solitons are weaker than in figure [21 even in cases with similar depletion: 
iV nc /iV tot ~ 0.2 in (d)-(f) as in figure |a)-(c). 

constant and N tot — > oo [EE]). Surprisingly, in the limit that g/N tot — > 0, 7 converges to 
a more negative value than that found in classical GPE simulations [figure Hfb)]. The 
classical soliton exhibits damping due to interaction with sound waves that originate 
from the imprinting process. Such interactions may represent dynamically unstable 
processes that require very weak numerical noise to be triggered which may be absent 
in GPE. Alternatively, the soliton trajectory may represent a state with non-classical 
statistics that does not converge to the classical GPE value, similarly to the atom 
number-squeezed states. Our simulations seem to indicate that the first case is a more 
likely explanation, since, e.g., an uncorrelated noise at each spatial grid point, with 
the magnitude ~ 10~ 4 weaker than quantum vacuum noise, reproduces the g/N tot — > 
TWA limit. 

4-2. Soliton core 

Many studies [521 [261 EB 1211 [29] suggest a link between dissipative soliton dynamics and 
filling of a soliton core. To investigate this, we calculate the ratio p s j pb of the minimum 
density in the soliton core to the background density around the soliton. For the classical 
soliton solution these variables are related by ([3]). We stress the distinction between the 
quantum expectation value of the soliton depth, which we infer from the expectation 
value (ps/Pb) found by locating the soliton and evaluating its depth in each stochastic 
realization, and the value of the density notch in the quantum expectation value (p(x)), 
obtained by ensemble averaging over densities in many realizations. We first calculate 
the quantum expectation value (p s /pb) at time t = 6.7 x 10~ 3 /u - very soon after 
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Figure 4. Quantum expectation value for atom density in soliton core divided by the 
background density (p s /pt) at T = 0. (a) at time t = 6.7 x 10~ 3 /w (immediately after 
imprinting). Larger quantum fluctuations indicate a shallower soliton core; (b) when 
the soliton crosses the point x — —0.4/ for the first time (lower/red line) and for the 
second time (upper/black line) . Larger quantum fluctuations now typically indicate a 
deeper soliton core that exhibits some filling as a function of time. 



the phase imprinting when the density dip is still forming and all solitons in different 
stochastic realizations are still approximately centered at x = 0. We approximate 
Pb by p(0.39Z) in each trajectory. Figure Ufa) shows how the filling of the emerging 
notch increases due to quantum fluctuations at T = 0, as density fluctuations obscure 
the soliton. However, by time t = 6 x 10~ 2 /u; after the formation of the density 
dip, the filling behaviour reverses and (p s /pb) becomes smaller in systems with large 
depletion. This increase in soliton depth due to quantum fluctuations is consistent with 
the aforementioned slowing down of the soliton velocities found in our previous study 
[2~T] in which case the increased phase uncertainty reduced the quantum expectation 
value of soliton speed (| cos(0/2)|) due to the nonlinear dependence of the speed upon 
the phase jump. The classical formula between the phase jump and soliton core density 
(J3]) also exhibits a negative curvature and we expect strong quantum fluctuations to 
reduce the quantum expectation value of soliton core density, reflected by (cos 2 (0/2)). 
This deepening of solitons is surprising considering recent studies [52], [26j [30] which 
suggest that a soliton core in individual runs may develop density peaks due to ground- 
state depletion, and we would naively expect shallower, faster solitons. Instead we find 
that the depth of a soliton in individual runs is dominated by large phase fluctuations 
which give rise to the opposite effect, so that, after the formation, a soliton core is deeper 
in systems with large depleted fraction. 

We now consider the behaviour of the soliton depth during the evolution in a 
harmonic trap. Figure IU(b) shows the quantum expectation value (p s /pb) as the soliton 
passes the point x = —0.4/ for the first and second time during its oscillation, which 
occur at times t ~ 4.8 and t ~ 8. In this case we may determine p s in every trajectory 
by the density at position x at the time when the soliton passes, and pb by averaging 
the density at the point x between consecutive passes of the soliton. We find that for 
individual trajectories, the soliton depth decreases for each pass; i.e., there is a slight 
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filling of the soliton core with time. This is accompanied by a gradual speeding of 
the solitons due to dissipation (7 < 0). Quantum fluctuations enhance filling of the 
soliton core. After formation, the soliton in systems with strong quantum fluctuations, 
however, is initially typically deeper than in the case of weak fluctuations, and figure 0] 
shows that they remain deeper despite experiencing greater filling. 

Numerically tracking the location of the phase kink in each individual realization, 
we have established that the expectation value for the depth of the soliton core is larger 
for larger values of N nc /N tot . Solitons may drift around due to quantum fluctuations 
and in different runs they are generally not located in the same spatial position at 
the same time - soliton position uncertainties grow as a function of time due to both 
initial velocity uncertainties and interactions of solitons with sound waves [21] . In an 
ensemble averaged density over many single realizations, corresponding to the quantum 
expectation value for the atom density (p(x,t)), soliton core density profiles appear 
broader and shallower due to different soliton locations in individual runs, as displayed 
in figure [5j The soliton position has a particularly large uncertainty in the presence of 
strong fluctuations and we find that, despite the fact that the expectation value of the 
soliton core is deeper in individual realizations with stronger quantum fluctuations, the 
density notches become shallower in (p(x,t)) when the fluctuations are enhanced. The 
flattening of (p(x,t)) due to the wandering solitons is similar to the previous results 
obtained using a Bogoliubov analysis [281 129] around a stationary dark soliton state 
showing that (p(x,t)) is smeared out due to fluctuations in the soliton positions. 

Matrix product state simulation of a dark soliton in a bosonic atomic gas in a lattice 
in the tight-binding approximation with unit filling was studied in references (3TJ [32] . 
From the expectation value of the atom density it was shown that an instantaneously 
imprinted phase kink with a vanishing density at the soliton core at the center of the 
lattice decays due to quantum fluctuations, as the soliton core gets filled with atoms. 
This is consistent with the aforementioned Bogoliubov analysis [28] 129] . Moreover, 
based on the time-evolution of the atomic pair-correlation function it was argued that 
the soliton core also in single-shot experimental realizations is filled. In the simulations 
the initial state of g 2 (0,k) = (alala ak) = 0, since the central site j = was empty 
(the location of the density dip of the soliton) and the other sites had approximately 
one atom. (Here denotes the annihilation operator for the atoms at the site k.) At 
later times g 2 became a flatter (and non-vanishing) function of k. Whether the filling 
of a soliton core in single-shot experiments in such a system duly happens, however, 
is inconclusive. This is because the evolution of g 2 can in many cases equally well 
be explained by a soliton whose core does not get filled in individual single-shot runs, 
but which randomly drifts along the lattice with the standard deviation of the soliton 
position increasing as a function of time. As a simple, idealized example, consider a 
lattice system where the long-range correlations may be ignored so that #2(0, k) ~ n nk- 
Here represents the atom density at the site k that is obtained by an ensemble average 
over many single-shot runs. The soliton is initially assumed to be located at the central 
site, with n = and rij = 1, for j 7^ 0, and g 2 (0, k) = for all k. At some later time 
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Figure 5. (a) Quantum expectation value for the soliton position (x(t)) (lines) 
and its uncertainty 8x (shaded area) for solitons in a harmonic trap at T = for 
(with increasing oscillation amplitudes) iVtot = 50 (magenta), iVtot = 440 (cyan) and 
AT tot = 8000 (negligible width); (b) as in (a), but with N to t = 900, T = 22fiuj/k B 
(magenta) and T — (cyan); (c) the soliton position uncertainty 5x at two hxed times 
ti ~ 1/uj (lower/black line) and £2 — 3/w (upper/red line) at T = 0, with iVtot varying 
between 40 and 8000; (d) as in (c) but varying T between and lQUw/ks for fixed 
Ntot = 900; (e)-(f) expectation value of atom density (p) at time t ~ 3.7/u;, obtained 
from I^vkI 2 by converting from the Wigner to the normally-ordered representation; 
(e) at T = curves correspond to iVtot = 50 (intermediate/red line), iVtot = 440 
(light/green line) and iV to t = 8000 (dark/black line); (f) for fixed N to t = 900 with the 
curves corresponding to T = (dark/black line) and T = IQUlo/Ub (light/red line). 
Increasing position uncertainty at low atom numbers create a filling effect when the 
density is averaged over many realizations. 

we may have nk = (L — 1)/L and #2(0, k) ~ (L — 1) 2 /L 2 , for all k, where L denotes 
the number of sites. This can represent an outcome where every single-shot run yields 
a constant atom density (L — 1)/L for all k, but equally well a situation where each 
individual run has a single soliton at a random location along the lattice, so that the 
probability of finding a soliton (with a vanishing density) at an arbitrary site k in each 
realization is 1/L and the probability of not finding a soliton at the same site (L — l)/L. 

Our numerical simulations of quantum dynamics of a soliton with the corresponding 
classical value of the imprinted phase jump C = ir in a lattice demonstrate a similar 
phenomenon (as shown in Section [6]): Individual stochastic realizations exhibit drifting 
soliton dynamics along the lattice due to quantum fluctuations. The ensemble average 
of the atom density and the pair-correlation function #2 become flatter as a function 
of time when the solitons have more time to drift over longer distances in the lattice. 
At the same time, however, the soliton cores in single-shot runs show very little, if any, 
effects of filling. 
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Figure 6. Individual stochastic realizations of \ijjw(%, t)\ 2 for collisions of two solitons 
at T — (top row), and corresponding trajectories of single solitons (bottom row) 
for fixed gN to t = WOhul, with i\T to t = 100 (left and centre) and 2V to t = 900 (right). 
At low atom number, for the case of strong quantum fluctuations, the solitons split 
after the collision. Comparison between the top and bottom rows shows how quantum 
fluctuations turn collisions inelastic, resulting in splitting solitons. 

5. Colliding solitons 

Recent studies of quantum dynamics of solitons in an optical lattice in the tight-binding 
limit with one mode function per lattice site have shown that soliton density may spread 
after colliding due to quantum fluctuations [31] [32] . We examine soliton collisions in 
the absence of a lattice, and consider the relationship of such inelastic collisions with 
the soliton structure in the presence of strong quantum fluctuations. 

We simulate the generation of two solitons using an optical imprinting potential (j3J) 
with two steps = 4166. 7hu at Xf, = —x a = 0.49/ which in a classical system would 
imprint phase kinks of C = 2 and C = 2tt — 2 with opposite velocities ±cos(0 c /2). 
We study the effects of quantum fluctuations on the imprinting of the phase kinks and 
the subsequent collisions of two such kinks (at about t ~ 0.2/u). As before, we do this 
by varying the effective interaction strength, corresponding to different values of the 
ground-state depletion, by keeping the nonlinearity N tot g fixed but adjusting the ratio 
g/N tot . Figure [6] shows l^jy^)! 2 f° r individual realizations of soliton collisions, and, for 
comparison, trajectories with the same stochastic realization, but using an imprinting 
potential that imprints only one of the solitons. We observe splitting/ spreading of the 
solitons as they emerge from collisions in the presence of large quantum fluctuations. 
Since the corresponding trajectories of single solitons do not exhibit this behaviour, 
we conclude that the splitting is induced by quantum fluctuations during the soliton 
collisions. Despite the spreading of the region of low atom density around the soliton 
core in individual realizations, the minimum density in the soliton core does not increase 
in these trajectories. We infer that one soliton splits into many solitons of comparable 
depth and velocity. 
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Figure 7. Individual stochastic realizations of (a^, *) | 2 for repeated collisions 
between different solitons. (a)-(c) Two solitons with (a) iVtat = 100 (b) TVtot = 1700 
and (c) Mot = 8000. (d)-(f) Six solitons with (d) and (e) 7Y tot = 100, (f) Mot = 900. 
Trajectories (d)-(e) show dramatic perturbations of the soliton motion, causing some 
solitons in (e) to disappear. 

Due to the dipolar motion in a harmonic trap, two solitons can experience repeated 
collisions. Recent experiment [43] showed seven repeated collisions of two solitons 
generated by merging two BECs. The parameters were close to those in our system: 
u ~ 333Hz, u r ~ 5592Hz and N tot ~ 1700 (gN tot ~ 205hul); or u ~ 364Hz, 
u r ~ 2563Hz and iVtot — 950 (gN tot ~ 50.1hul). The experiment was performed in an 
elongated ID trap, but without sufficiently strong transverse confinement that would 
have suppressed radial density fluctuations, and quantum fluctuations are not expected 
to be important. In our ID system we find little quantum effects at similar atom 
numbers (other than a reduction in the soliton velocities). However, after reducing the 
atom number to iVtot = 100, we observe soliton trajectories becoming perturbed [figure 
UJa)] or disappearing after several collisions that may be detectable in an experiment 
using a tighter transverse confinement than the one in |43j . 

Increasing the number of solitons from two introduces the possibility of far richer 
dynamics. Systems of three or more harmonically trapped bright solitons may exhibit 
chaotic dynamics; the dynamics are not integrable due to the interplay between 
the harmonic motion and the soliton interaction [53]. The similarity between the 
interactions of bright solitons and those of dark solitons suggests that multiple dark 
solitons in a harmonic potential may also display chaotic dynamics. Multiple soliton 
collisions may also appear, e.g., in self- focusing and revival dynamics of BECs [54J. 
As a signature of quantum fluctuations, we find that in such systems, solitons are 
more susceptible to disappearing. We use the imprinting potential (j4|) to imprint six 
phase kinks all with different initial positions and velocities. Figure [7(d)- (f) shows 
|-0w(^)| 2 for individual realizations. The enhanced dissipative effects visible in[7Jd)-(e) 
where solitons trajectories become rapidly perturbed or disappear, may be related to 
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dynamical instability of the system associated with chaotic dynamics, or merely due to 
the increased number of collisions experienced by each soliton in a short period of time. 



6. Soliton in a combined harmonic trap and optical lattice 

The classical nonlinear dynamics of a phase kink in a combined harmonic trap and 
optical lattice exhibits dynamical instabilities and the soliton may dissipate energy via 
sound emission even without quantum and thermal fluctuations [55} 56J. The interplay 
between quantum fluctuations and the dynamical instabilities of the corresponding 
classical system was previously studied by us [21]. We found that the fluctuations can 
enhance the effect of instabilities, in some situations resulting in a surprising reduction 
in the position uncertainty of the solitons. Single-shot soliton trajectories displayed 
splitting and disappearing solitons in an optical lattice, similar to those following 
collisions, studied in section |5j Here we will focus on solitons that retain their integrity 
during the evolution. We consider very slow solitons that would have zero initial velocity 
in the classical case. Such a system exhibits strong fluctuation-induced effects and we 
calculate quantum statistical properties of the soliton dynamics, such as atom number 
fluctuations in an individual sites, waiting time distributions for the soliton to exit the 
initial site, atom populations and pair-correlations. 

In the simulations of imprinting and dynamics we include in ([5]) the optical lattice 
potential V exp = mu 2 x 2 /2 + sE r sin 2 (7nc/(i), where E r = h 2 n 2 /2md 2 is the lattice photon 
recoil energy and d is the lattice period. Here we set d = ttI/A and s — 1. After the 
simulations, in order to translate Wigner distributed stochastic variables to normally 
ordered expectation values, following the approach in [91 [TT] we define the ground-state 
operators Oj for the individual lattice sites as 

aj = / dxip*(x)il} W (x,t) (14) 

J jth well 

where ipj is the Gaussian ground-state wavefunction (Wannier function) of the jth well 
and ipw{ x , t) is the numerically integrated full Wigner wavefunction. The population of 
the jth site is thus (aja-,) = (a*aj)w — 1/2, and the number fluctuations in the jth site 



5 nj 



({a]aj) 2 ) - (a]ajf 



1/2 



((a* aj ) 2 ) w - {a*aj) 2 w - - 



1/2 



(15) 



The pair-correlation function between the fcth and the Oth site (the central site in which 
the soliton is imprinted) reads 

5f 2 (0, k) = (alala k a ) (16) 

— ( a o a l a k a o)w — - ^(«o a o)vK + {o* k a k ) w — — 5 k (a* k ak)w — ^ • 

We apply an imprinting potential that would, in the classical (GPE) case, prepare 
a phase kink of C = 7r (a soliton with zero velocity). This state is unstable [55], as small 
oscillations in a lattice site become amplified by sound emission so the soliton can escape 
from the central site and begin to drift along the lattice [56]. Quantum fluctuations seed 
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Figure 8. Solitons in a combined harmonic trap and optical lattice, with iVtot = 900 
and T = 0. (a) The quantum expectation value of the time r s that it takes for a soliton 
to exit the central lattice site and its quantum uncertainty (error bar) vs the phase 
jump (p c ; (b) quantum expectation values for the soliton position (x(t)) (lines) and its 
uncertainty 5x (shaded areas) for <f> c = tt (magenta) and 4> c — 2 (yellow); (c)-(d) the 
Wigner density \ipw(x, t)\ 2 for individual realizations of soliton dynamics with (fi c = 7r; 
(e) the correlation function 32(0, x); (f) the density expectation value in the lowest 
band, (e) and (f) both show distributions becoming fiat due to randomly drifting 
solitons along the lattice without solitons in individual realizations disappearing. 

this soliton motion as demonstrated in figure E^e)-(f) where the solitons exit the central 
cite at different times and in different directions. Figure [S^a) shows that the quantum 
uncertainty in the time t s that it takes for a soliton to exit the central site becomes 
comparable to the quantum expectation value r s . 

The population expectation value for the central lattice site (no) initially shows a 
soliton oscillating at the frequency of the lattice site, indicating overlapping oscillations 
between this site and the adjacent sites. This is followed by a rapid increase in central 
site occupation (no) after time t ~ 4/u when there is a high probability that the soliton 
has left the site. This behaviour is reflected in the expectation value of the atom density 
in the lowest energy band at x = [figure 1ST)]. We also find an initial oscillatory 
behaviour followed by a rapid increase in the pair correlation function [figure El^e)]. 
Previous studies of solitons in lattices with small atomic occupations J3TJ [32] cited such 
an increase as evidence of filling of a soliton core in individual realizations, but here in 
large atom number systems we find a similar effect caused by solitons drifting along the 
lattice, while a soliton core in none of the individual realizations is significantly filled. 

As well as affecting the expectation values (no) and g2, the drifting behaviour 
of solitons has dramatic effects on the atom number fluctuations. Figure shows 

1 /2 

the number squeezing parameter in the central lattice site 5no/n corresponding to 
different <ft c generated by changing the imprinting time r in the potential (jlj). When the 
fluctuations obey a Poisson distribution, 5no/n = 1. The atom number fluctuations 
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Figure 9. Atom number squeezing parameter in the central lattice site for a soliton in 
a combined harmonic trap and optical lattice, with iVtot — 900 and T = 0. The phase 
notch across the soliton is imprinted by a potential which in the corresponding classical 
system imprints the phases (a) </> c = 1.3 (fast soliton), (b) <f> c = 2.6, (c) <f> c = 3.0 (d) 
4> c = 7r (classical soliton has zero velocity). Soliton crossing x = is identified by 
super-Poisson atom number fluctuations. At lower soliton velocities, the sharp peaks 
are broadened due to the quantum uncertainty in the soliton crossing time. 



of a ground-state BEC in an optical lattice are sub-Poissonian (squeezed), such that 
8no/n < 1 [9| [TTj . For fast solitons, the soliton position uncertainty is small [figure [8] 
(b)], soliton cores between different realizations overlap, so there are predictable times 
when the soliton at high probability is not in the central lattice site. During these times 
number statistics in the site remain squeezed. The crossing of the soliton at x = 
can be recognized as a sharp super-Poissonian peak in the dynamics of the squeezing 
parameter. For slower solitons, the position uncertainty of the soliton is so large that 
there are large atom number fluctuations in the central lattice site due to fluctuations 
in the crossing times of the solitons at x — 0. 
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